//
// Created by Ivor_ on 2023/2/8.
//

#include "logisticLike.h"

static double var1;
static double var2;
static double var3;
static double var4;

static double linearFunction2PF(double xx);

double* logistic(int iterDim, double* iterVar, int argCount, double* args) {
    double* logisticIterVars = (double*)calloc(iterDim, sizeof(double));
    for (int i = 0; i < iterDim; ++i) {
        logisticIterVars[i] = args[i] * (iterVar[i]-pow(iterVar[i], 2.));
    }
    return logisticIterVars;
}

double* partialDerivativeLogistic(int partialVarIndex, int iterDim, double* iterVar, int argCount, double* args) {
    double* logisticIterDerVars = (double*)calloc(iterDim, sizeof(double));
    for (int i = 0; i < iterDim; ++i) {
        logisticIterDerVars[i] = args[i] * (iterVar[i]-pow(iterVar[i], 2.));
    }
    logisticIterDerVars[partialVarIndex] = (1. - 2.*iterVar[partialVarIndex]) * args[partialVarIndex];
    return logisticIterDerVars;
}

void plotIterChanges(char* fileName, int plotDim, int plotIterMinIndex, int plotIterMaxIndex, double yMin, double yMax, double*(*iterFun)(), int iterDim, double* iterVar, int argCount, double* args) {
    PIXEL whitePixel = {.red = 0xff, .green = 0xff, .blue = 0xff};
    PIXEL blackPixel = {.red = 0x0, .green = 0x0, .blue = 0x0};
    STDbmp* picture = bmpBO().new(fileName, 1920u, 1080u, whitePixel);
    SCALE scale = {.xStart = plotIterMinIndex - .2*(plotIterMaxIndex - plotIterMinIndex), .xEnd = plotIterMaxIndex + .2*(plotIterMaxIndex - plotIterMinIndex), .yStart = yMin, .yEnd = yMax};
    for (int i = 0; i < plotIterMaxIndex; ++i) {
        if (i >= plotIterMinIndex) {
            bmpPlt().ptsIter(picture, scale, 7, 7, blackPixel, 1, 2, 1.*(i - plotIterMinIndex), iterVar[plotDim]);
            if (i == plotIterMinIndex) {
                var4 = 1.*(i - plotIterMinIndex);
                var3 = iterVar[plotDim];
            } else {
                var2 = var4;
                var1 = var3;
                var4 = 1.*(i - plotIterMinIndex);
                var3 = iterVar[plotDim];
                bmpPlt().invFun(picture, scale, linearFunction2PF, var1, var3, 3, 3, blackPixel);
            }
        }
        iterVar = iterFun(iterDim, iterVar, argCount, args);
    }
    bmpBO().close(picture);
}

double feigenbaum(char* fileName, int plotDim, int objArgIndex, double objArgMin, double objArgMax, double yMin, double yMax, int iterN1, int iterN2, double*(*iterFun)(), int iterDim, double* iterVar, int argCount, double* args) {
    PIXEL whitePixel = {.red = 0xff, .green = 0xff, .blue = 0xff};
    PIXEL blackPixel = {.red = 0x0, .green = 0x0, .blue = 0x0};
    STDbmp* picture = bmpBO().new(fileName, 1920u, 1080u, whitePixel);
    SCALE scale = {.xStart = objArgMin - .2*(objArgMax - objArgMin), .xEnd = objArgMax + .2*(objArgMax - objArgMin), .yStart = yMin, .yEnd = yMax};
    double* bifurcationData = (double*)calloc(6, sizeof(double));
    bifurcationData[0] = 1.;
    bifurcationData[1] = objArgMin;
    bifurcationData[2] = 1.;
    bifurcationData[3] = objArgMin;
    bifurcationData[4] = 1.;
    bifurcationData[5] = objArgMin;
    for (int j = 0; j < (unsigned int)((objArgMax-objArgMin) / (scale.xEnd-scale.xStart) * 1920.); ++j) {
        args[objArgIndex] = objArgMin + j*(scale.xEnd-scale.xStart)/1920.;
        double temp, counter;
        for (int i = 0; i < iterN2; ++i) {
            if (i == iterN1) {
                temp = iterVar[plotDim];
                counter = 0.;
            }
            iterVar = iterFun(iterDim, iterVar, argCount, args);
            counter += 1.;
            if (i >= iterN1) {
                bmpPlt().ptsIter(picture, scale, 3, 3, blackPixel, 0, 2, objArgMin + j * (scale.xEnd - scale.xStart) / 1920., iterVar[plotDim]);
                if (fabs(bifurcationData[5] - args[objArgIndex]) > (scale.xEnd-scale.xStart)/1920. && fabs(temp - iterVar[plotDim]) < 1.e-12*(yMax - yMin) && bifurcationData[4] < (iterN2 - iterN1)/2) {
                    if (counter > bifurcationData[4]) {
                        bifurcationData[0] = bifurcationData[2];
                        bifurcationData[1] = bifurcationData[3];
                        bifurcationData[2] = bifurcationData[4];
                        bifurcationData[3] = bifurcationData[5];
                        bifurcationData[4] = counter;
                        bifurcationData[5] = args[objArgIndex];
                    }
                    counter = 0.;
                }
            }
        }
    }
    bmpBO().close(picture);
    double feigenbaumDelta = (bifurcationData[3]-bifurcationData[1]) / (bifurcationData[5]-bifurcationData[3]);
    free(bifurcationData);
    return feigenbaumDelta;
}

void lyapunov(char* fileName, int plotDim, int objArgIndex, double objArgMin, double objArgMax, double yMin, double yMax, int iterNum, double*(*iterFun)(), double*(*iterDerFun)(), int partialIndex, int iterDim, double* iterVar, int argCount, double* args) {
    PIXEL whitePixel = {.red = 0xff, .green = 0xff, .blue = 0xff};
    PIXEL blackPixel = {.red = 0x0, .green = 0x0, .blue = 0x0};
    STDbmp* picture = bmpBO().new(fileName, 1920u, 1080u, whitePixel);
    SCALE scale = {.xStart = objArgMin - .2*(objArgMax - objArgMin), .xEnd = objArgMax + .2*(objArgMax - objArgMin), .yStart = yMin, .yEnd = yMax};
    bmpBO().drawLine(picture, (unsigned int)(-1. * picture->height * yMin / (yMax-yMin)), 7, picture->width/10, blackPixel);
    double x0, y0;
    for (int i = 0; i < (unsigned int)((objArgMax-objArgMin) / (scale.xEnd-scale.xStart) * 1920.); ++i) {
        args[objArgIndex] = objArgMin + i*(scale.xEnd-scale.xStart)/1920.;
        double yy = 0.;
        for (int j = 0; j < iterNum; ++j) {
            iterVar = iterFun(iterDim, iterVar, argCount, args);
            yy += log(fabs(partialDerivativeLogistic(plotDim, iterDim, iterVar, argCount, args)[plotDim]))/iterNum;
        }
        if (i == 0) {
            x0 = objArgMin + i*(scale.xEnd - scale.xStart)/1920.;
            y0 = yy;
        }
        bmpPlt().ptsIter(picture, scale, 3, 3, blackPixel, 1, 4, x0, y0, objArgMin + i*(scale.xEnd - scale.xStart)/1920., yy);
        x0 = objArgMin + i*(scale.xEnd - scale.xStart)/1920.;
        y0 = yy;
    }
}

static double linearFunction2PF(double xx) {
    return (xx - var3)/(var1 - var3)*var2 + (xx - var1)/(var3 - var1)*var4;
}
